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Abstract 

A combination of classical molecular dynamics (MD) and ab initio Car-Parrinello molecular 
dynamics (CPMD) simulations is used to investigate the adsorption of water on a free amorphous 
silica surface. From the classical MD SiC>2 configurations with a free surface are generated which 
are then used as starting configurations for the CPMD. We study the reaction of a water molecule 
with a two-membered ring at the temperature T = 300 K. We show that the result of this reaction 
is the formation of two silanol groups on the surface. The activation energy of the reaction is 
estimated and it is shown that the reaction is exothermic. 
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1 



I. INTRODUCTION 



The interaction of water with surfaces of amorphous silica is of great technologi- 
cal interest 111 M, and thus numerous studies have been devoted to this issue ranging 

to numerical methods usin g bo th ab initio tech- 

17, 18, Q. In these 



Ft spectroscopy 3, 

flfiQQE, 



from IR spectrosco 
niques 



14L Il5j and classical MD simulations 
studies, the possible existence of small-membered Si-0 rings on Si02 surfaces has been in- 
tensely discussed since such rings are expected to be the reactive centres for the interaction 
with water and other molecules. It is well-known that water may dissociate on SiC>2 surfaces 
resulting in the formation of silanol (Si-OH) groups. These silanol groups can be detected 
in spectroscopic experiments via the Si-OH stretching mode near 3750 cm" - 1 120. 21 



In particular it is believed that the silanol groups are a result of the interaction of water 
molecules with small-membered rings 

mm. 

In this paper we investigate the reaction of a water molecule with a two-membered ring 
on an amorphous silica surface. To this end, we combine classical MD simulations u sing 
the so-called BKS potential Q for silica and ab initio Car-Parrinello MD (CPMD) 24 1. 
This combination of methods has recently also been used for investigating silicates in the 



bulk 125 



261 ] and of free silica surfaces 27j . The idea in our case is to generate configurations 



with a free surface by "BKS-MD" which are then the starting configurations for CPMD runs. 
The advantages of this methodology are two-fold: On the one hand, in the classical MD 
amorphous systems can be relaxed on a nanosecond time scale whereas typical time scales 
for the CPMD are only of the order of several picoseconds. On the other hand, the CPMD 
provides a much more realistic modeling of atomistic systems. This is of special importance 
in the case of interactions of water with an Si02 surface since it is difficult to model these 
interactions realistically by classical potentials. 

The rest of the paper is organized as follows: In the next section we describe how we have 
prepared a free silica surface by classical MD and give the details of the CPMD simulations. 
Then, the results for the water interaction with the two-membered ring are presented in 
Sec. 3. Finally, we summarize the results. 
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II. PREPARATION OF THE SURFACE AND DETAILS OF THE SIMULATIONS 



In order to study the reaction of a silica surface with a water molecule, we first have 
to prepare SiC>2 configurations with a free surface. As we have already mentioned in the 
Introduction we use for this a combination of classical MD and CPMD. 



2311 which 



As a model potential for the classical MD we use the so-called BKS potential 
is a simple pair potential of the following form: 

cf>(r) = ^^+A aP exp(-B aP r)-^ a,/?e[Si,0], (1) 

where r is the distance between an atom of type a and an atom of type (3. The effective 
charges are go — —1-2 and gsi — 2.4, and the parameters A a p, B a p, and C a p can be found 
in the original publication. They were determined by using a combination of ab initio 
calculations and classical lattice dynamics simulations. The long-ranged Coulomb forces 
(and the potential) were evaluated by means of the Ewald summation technique. As an 
integrator for the simulation we used the velocity form of the Verlet algorithm with a 
time step of 1.2 fs. 

In order to investigate a system with free surfaces one could consider a film geometry, i.e. a 
system with periodic boundary conditions (PBC) in two directions and free boundaries in the 
third direction. Unfortunately, this is not a very good approach since the Ewald summation 
technique for the long ranged Coulomb interactions becomes inefficient in this case. This 
stems from the fact that the Fourier part of the Ewald sums can no longer be calculated by 
a single loop over the number of particles N as in the case of PBC in all three directions, 



29 



30, 311. Therefore, we have 



but one has to compute a double loop that scales with N 2 
adopted in this work a different strategy in that a system with PBC in all three directions 
was simulated containing an empty space Az in ^-direction. 

The preparation of the interface is done by the following steps: i) We start with a system 
at T = 3400 K with PBC in all three directions (box dimensions: L x — L y — 11.51 A and 
L' z = 23 A). This system is fully equilibrated within 1 ns (see also below). As a result, 
configurations are obtained that indicate a realistic modelling of Si02 with the empirical 
BKS potential, Eq. Q: A tetrahedral Si-0 network is formed containing defects (e.g. five- 
fold coordinated silicon atoms) which are expected at a temperature as high as T = 3400 K. 
However, no artificial bonds such as Si-Si or 0-0 bonds are found in the resulting network 
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structures, ii) We cut the system perpendicular to the ^-direction into two pieces. This 
is only done for oxygen-silicon bonds such that we get only free oxygen atoms at this 
interface, iii) These free oxygen atoms are now saturated by hydrogen atoms. The place 
of these hydrogen atoms is chosen such that each of the new oxygen-hydrogen bonds is 
in the same direction as the oxygen-silicon bonds which were cut and have a length of 
approximately 1 A. The interaction between the hydrogen and the oxygen atoms as well 
as the silicon atoms are described only by a Coulombic term. The value of the effective 
charge of the hydrogen atoms is set to 0.6 which ensures that the system is still (charge) 
neutral, iv) We permanently freeze atoms which have a distance from the interface that 
is less than 4.5 A, whereas atoms that have a larger distance can propagate subject to the 
force field, thus generating a mobile layer of 14.5 A. v) We add in z-direction an empty 
space Az = 6.0 A . With this sandwich geometry we can use periodic boundary conditions 
in all three directions. We have made checked that the value of Az is sufficiently large to 
ensure that the structure of the free surface is not affected by the immobilized part of the 
system j^. Finally, we have a system of 91 oxygen, 43 silicon, and 10 hydrogen atoms in a 
simulation box with L x = L y = 11.51 A and L z m 25 A. 

We have done 125 independent runs where we have equilibrated the system for 1 ns. For 
the present study we have selected from these runs one configuration that has exactly one 
two-membered ring on its surface. This configuration was then quenched to T — 300 K 
using a cooling rate of 2.8 x 10 12 K/s where it was further relaxed for 0.25 ns. Subsequently 
the final configuration of this run has been used as a starting point for the CPMD. 

In the CPMD, we used conventional pseudopotentials for silicon and oxygen and the 
BLYP exchange functions 33, 3^|. The electronic wave-functions were expanded in a plane 
wave basis set with an energy cutoff of 60 Ry and the equations of motion were integrated 
with a time step of 0.085 fs for 0.2 ps. In the analysis of the CPMD run only those config- 
urations were taken into account that were produced later than 5 fs after the start of the 
CPMD run in order to allow the system to equilibrate at least locally j^ . 

III. RESULTS 

On intermediate length scales the surface structure of SiOo can be well distinguished from 
its bulk structure by the ring distribution [35]. A ring is defined as the shortest closed loop 
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a) 



n=2 



b) 



n=3 




FIG. 1: Geometry of a) two-member ed rings and b) three-membered rings as obtained from BKS- 
MD (indicated by solid lines) and from CPMD (indicated by the dashed lines). Si atoms are shown 
as light circles and O atoms as small dark circles. For the sake of clarity of the figure, distances 
and angles are not to scale, and hence the difference between BKS-MD and CPMD in the figure 
is somewhat exaggerated in comparison with reality. 



of n consecutive Si-0 elements. The Si02 surfaces that were obtained from the BKS-MD at 
T = 3400 K showed a relatively high concentration of 2-membered rings, i.e. edge-sharing 
tetrahedra, in contrast to the bulk in which, even at high temperatures, such structures 
are basically absent. The relaxation of the surfaces by CPMD reduced the number of 2- 
membered rings significantly at T — 3400 K. At T = 300 K, the time scale on which the 
system is relaxed by means of CPMD (0.2 ps) is not sufficient to cause a rearrangement 
of ring structures and thus the number of 2-membered rings does not change during the 



CPMD 



3J. However, both at T = 3400 K and at T = 300 K, the CPMD yielded a 



slightly different geometry of 2-membered rings which is illustrated in Fig.^,. For BKS-MD 
as well as CPMD one finds a nearly planar geometry of two-membered rings. Whereas in 
CPMD a nearly quadratic shape with an O-Si-0 angle around 90° is seen, the BKS-MD 
leads to a trapezoid geometry with an O-Si-0 angle around 80°. This is due to the fact 
that the distance between nearest Si-0 and 0-0 neighbors is larger in CPMD while the 
distance between the silicon atoms in the 2-fold ring is approximately around 2.5 A in both 
methods. A similar behavior can be seen for the 3-fold rings (see Fig. [Qd) which exhibit 
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FIG. 2: Snapshot that illustrates the preparation of the water molecule on a two-fold ring. Note 
that Si(l), Si(2), 0(1), and 0(2) are part of the Si-0 network. 



also a nearly planar geometry. Also, in this case BKS-MD and CPMD show a similar Si-Si 
distance (between nearest Si neighbors) around 3.0 A while the Si-0 distance is slightly 
larger in the CPMD. As a consequence, in the BKS-MD, the typical O-Si-0 angle is 99° 
whereas in the CPMD an O-Si-0 of 107° is found which is close to the angle of 109.47° in 
an ideal tetrahedron. 

Small-membered rings are of particular interest if one considers the reaction of a water 



molecule with an Si02 surface. A 



possib 



e reaction mechanism which was put forward 



by molecular orbital calculations 0, [lfj, |ll( implies the disruption of a Si-O-Si structural 
element followed by the formation of two silanol groups (Si-O-H). Thereby a Si-0 bond 
has to be broken. For this reason, 2- or 3-membered rings are of particular interest since 
these structures have a high local internal stress and hence their bonds can be broken more 
easily than the one in larger rings. Therefore, we consider in the following the reaction of 
H 2 molecule with a 2-membered ring on a silica surface. To this end, we have selected a 
configuration from our BKS-CPMD simulations at 300 K that exhibit exactly one two-fold 
ring on its surface. 

The initial position that we have chosen for the H2O molecule on the 2-membered ring 
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FIG. 3: Total energy E'tot an d potential energy Spot during the interaction of the silica surface with 
the H2O molecule (note the logarithmic time axis). The vertical dashed lines at t± = 0.0045 ps, 
t2 = 0.009 ps, and t% = 0.045 ps mark the times at which the snapshots in Fig. 0] are shown. 

is illustrated in Fig. |21 The oxygen atom of the water molecule [0(3)] was placed at a 
distance of 1.52 A from one of the Si atoms [Si(l)] which is slightly smaller than the Si-0 
bond lengths in the ring of 1.64 A. The two hydrogen atoms H(l) and H(2) were placed 
at a distance of 1 A from 0(3) (and H(l) at a distance of 1 A from 0(1)) such that the 
H(l)-0(3)-H(2) angle corresponds to 109°. Further details can be extracted from Fig. |21 
We shall see in the following that the small 0(l)-0(3) distance (1.63 A compared to 2.2 A in 
bulk silica) drives the breaking of a Si-0 bond and the formation of two silanol groups. 

The initial condition for H2O on the two-membered rings is of course artificial and leads 
to enormous forces between the atoms. In order to avoid an explosion of the system we 
have introduced a thermostat. If the temperature becomes larger than 380 K or lower than 
220 K, the velocities of the particles are rescaled such that the apparent temperature of 
the system is 300 K. Figure El shows the potential and the total energy during the reaction 
whereby the time t = corresponds to the initial configuration. Obviously, the potential 
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a) b) c) 




FIG. 4: Snapshots of the reaction of the H2O molecule with the 2-fold ring at a) t = t%, b) t = t%, 
and c) t = t$ (ti, and £3 as indicated in Fig. |3J). 

energy decreases rapidly such that after about 0.1 ps a constant is reached. Snapshots of 
the reaction at t\ = 0.0045 ps, t 2 = 0.009 ps, and £3 = 0.045 ps are shown in Fig. 0] (note 
that ti, t 2 , and £3 are also marked as vertical dashed lines in Fig. EJ). Indeed, at time £3 the 
end configuration with two silanol groups is formed. The first step that leads to the silanol 
groups is the breaking of one of the H-0 bonds and of the enforced 0(3)-0(l) bond. This 
has happened at t = t\ whereby the potential energy is already quite close to the equilibrium 
value (see Fig. EJ). In a second step then a Si-0 bond is broken (see snapshot at t = t 2 ) 
followed by an increased separation of the two silicon atoms. 

Now we want to study in more detail the formation of two silanol groups as a result of 
the reaction of the water molecule with a two-fold ring. This reaction can be written as 
follows: 




+ 5E 



(2) 

In the following we will estimate the activation energy A and also 5E. A positive or negative 
sign of 5E indicates whether the reaction is exothermic or endothermic, respectively. We use 
again the starting configuration shown in Fig. El But now we completely freeze in the SiC>2 
system and we move only the water molecule at a temperature of 30 K. As the inset of Fig. |3] 
shows, the potential energy is monotonously decaying while the water molecule relaxes on 
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FIG. 5: Change in the distance Sr between the indicated atoms (that are defined in the two 
snapshots at times t\ = 0.044 ps and t\ = 0.047 ps) during the simulation at T = 30 K where only 
the H2O molecule was moved whereas the other atoms were fixed. The inset shows the potential 
energy for the latter simulation (note that E po t is shifted by 1640 a.u.). 

the surface. After about 0.1 ps a constant value of E pot is reached where the water molecule 
has evaporated from the surface. Also shown in Fig. |5] are the distances Sr of the water 
atoms from the initial position as a function of time. We see that, slightly before the H 2 
evaporates, the 0(3)-0(l) distance and the 0(3)— Si(l) distance are increasing while the 
0(3)-H(l) distance is decreasing to the ideal O-H distance in a water molecule around 
1 A. Of importance for the following, are the configurations at the times t\ = 0.044 ps 
and t\ = 0.047 ps (marked as vertical dashed lines in Fig. |5J) that we will denote by SCI 
and SC2, respectively. Using SCI as a starting configuration for a CPMD run, where the 
whole system is again simulated at T = 300 K, leads to the formation of two silanol groups, 
whereas the use of SC2 yields the evaporation of the water molecule from the surface. As 
we can infer from the snapshots in Fig. the main difference between SCI and SC2 is a 
slightly different 0(3)-H(l) distance. 

In Fig. El we plot the potential energy for the two CPMD runs starting from the configu- 
rations SCI and SC2. Since both runs yield initially about the same value for E pot we can 
extract both A and SE from the plot. First, from the difference of the saturated values for 
t > 0.1 ps (marked as horizontal lines in Fig.EJ) we obtain SE = i?(Si0 2 +H 2 0) — £'(2SiOH) « 



9 




-4 -1 -? -1 

10 10 10 10 

t[ps] 

FIG. 6: Potential energy as a function of time for the simulations at T = 300 K using the start 
configuration SC2 (dashed line) and SCI (solid line). The snapshots correspond to the end config- 
urations of the two runs as indicated. In both figures the two hydrogen atoms are shown as small 
white spheres at the bottom of the simulation box. 

1.6 eV which indicates that the reaction is exothermic. The activation energy A can be es- 
timated from the difference of the initial and the final value of the SC2 run where the water 
molecule evaporates from the surface. For the latter process an HO H bond and a Si-0 
bond have to be broken and thus the binding energy of both bonds contribute to A. From 
Fig. El the value A 0.9 eV can be extracted. Whereas the distances of the 0-0 and Si-0 
neighbors in the two-membered ring are significantly smaller than those found in the bulk 
(see above), the same bonds are very close to the bulk values in the final structure with the 
two silanol rings. Moreover, the O-Si-0 angle in the silanol groups has a value around 108°, 
i.e. close to the value found in an ideal Si04 tetrahedron. These findings may explain the 
observation of an exothermic reaction. 

In a recent publication, Masini and Bernasconi also studied the adsorption reaction 
of a water molecule with a two-membered ring on an amorphous silica surface by CPMD. 
They find 5E ~ 1.7 eV, i.e. a value which is very similar to the one found in this work. 
However, they did a constrained MD simulation to estimate the activation energy. As a 
reaction coordinate they used either the distance between a silicon atom in the ring and the 
oxygen atom of the water molecule (called path A in Ref. Q|) or the distance between an 
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oxygen atom of the ring and a hydrogen atom of the water molecule (called path B). As a 
result A = 1.1 eV is obtained for path A and A = 0.32 eV for path B. In this work, we find 

n 

an activation energy which is closer to path A in Ref. [12J. 



IV. SUMMARY 

Using a combination of classical MD and CPMD we have investigated the reaction of 
water with a free amorphous silica surface. We have seen that the reaction of a water 
molecule with a two-membered ring, leading to the formation of two silanol groups on the 
SiC>2 surface, is an exothermic reaction with an activation energy around 0.9 eV. That the 
reaction is exothermic is in agreement with other recent numerical studies, and may be 
explained by the change of bond lengths and angles towards values found in bulk SiC>2 while 
the silanol groups form. 

Acknowledgements 

We are grateful to Michele Parrinello and Gloria Tabacchi for the introduction of one of us 
(C. M.) into CPMD. Without their support the present work would not have been possible. 
We acknowledge financial support by the SCHOTT Glaswerke Fond, the DFG under SFB 
262, the BMBF under grant No. 03N6015, and the European Community's Human Poten- 
tial Program under contract HPRN-CT-2002-00307, DYGLAGEMEM. J. H. acknowledges 
financial support from the DFG under grants HO 2231/2-1/2. We thank the NIC Julich for 
a generous grant of computing time. 



[1] A. P. Legrand, The Surface Properties of Silica, Wiley, New York, 1998. 

[2] R. K. Her, The Chemistry of Silica, Wiley, New York, 1979. 

[3] B. A. Morrow, I. A. Cody, J. Phys. Chem., 1976, 80, 1995. 

[4] T. A. Michalske, B. C. Bunker, J. Appl. Phys., 1984, 56, 2686. 

[5] B. C. Bunker, D. M. Haaland, K. J. Ward, T. A. Michalske, W. L. Smith, J. S. Binkley, C. 

F. Melius, C. A. Balfe, Surf. ScL, 1989, 210, 406. 

[6] L. H. Dubois, B. R. Zegarski, J. Am. Chem. Soc, 1993, 115, 1190. 



11 



[7] L. H. Dubois, B. R. Zegarski, J. Phys. Chem., 1993, 97, 1665. 

[8] A. Grabbe, T. A. Michalske, W. L. Smith, J. Phys. Chem., 1995, 99, 4648. 

[9] M. O'Keeffe, G. V. Gibbs, J. Chem. Phys., 1984, 81, 876. 

[10] T. Uchino, Y. Tokuda, T. Yoko, Phys. Rev. B, 1998, 58, 5322. 

[11] I. S. Chuang, G. E. Maciel, J. Am. Chem. Soc, 1996, 118, 401. 

[12] P. Masini, M. Bernasconi, J. Phys.: Condens. Matter, 2002, 14, 4133. 

[13] M. H. Du, L. L. Wang, A. Kolchin, H. P. Cheng, Eur. Phys. J. D, 20 03 , 84, 323. 

[14] M. H. Du, A. Kolchin, H. P. Cheng, J. Chem. Phys., 2003, 119, 6418. 

[15] M. H. Du, A. Kolchin, H. P. Cheng, J. Chem. Phys., 2004, 120, 1044. 

[16] E. B. Webb, S. H. Garofalini, J. Non-Cryst. Sol, 1998, 226, 47. 

[17] V. A. Bakaev, W. A. Steele, J. Chem. Phys., 1999, 111, 9803. 

[18] M. Wilson, T. R. Walsh, J. Chem. Phys., 2000, 113, 9180. 

[19] T. R. Walsh, M. Wilson, A. P. Sutton, J. Chem. Phys., 2000, 113, 9191. 

[20] M. L. Hair, Infrared Spectroscopy in Surface Chemistry, Dekker, New York, 1967. 

[21] P. R. Ryason, B. G. Russell, J. Phys. Chem., 1975 , 79, 1276. 

[22] B. A. Morrow, A. J. McFarlan, J. Non-Cryst. Sol, 1990, 120, 61. 

[23] B. W. H. van Beest, G. J. Kramer, R. A. van Santen, Phys. Rev. Lett. 1990, 64, 1955. 

[24] R. Car, M. Parrinello, Phys. Rev. Lett. 1985, 55, 2471. 

[25] M. Benoit, S. Ispas, P. Jund, R. Jullien, Eur. Phys. J. B 2000, 13, 631. 

[26] S. Ispas, M. Benoit, P. Jund, R. Jullien, Phys. Rev. B, 2001, 64, 214206. 

[27] C. Mischler, W. Kob, K. Binder, Comp. Phys. Comm., 2002, 147, 222. 

[28] M. P. Allen, D. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. 

[29] D. E. Parry, Surf. Sci., 1975, 49, 433. 

[30] D. E. Parry, Surf. Sci, 1976, 54, 195. 

[31] S. W. de Leeuw and J. W. Perram, Physica A, 1982, USA, 546. 

[32] C. Mischler, Molekulardynamik-Simulation zur Struktur von Si02-Oberflachen mit adsor- 

biertem Wasser, Ph.D. Thesis, Mainz, 2002. 

[33] N. Trouiller, J. L. Martins, Phys. Rev. B, 1991, 43, 1993. 

[34] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B, 1988, 37, 785. 

[35] A. Roder, W. Kob, K. Binder, J. Chem. Phys., 2001, 114, 7602. 



12 



